#include <deque>
#include <map>
#include <vector>
#include <algorithm>

#include "../sdf/dynarray.h"
#include "../contree/disjointset.h"
#include "dctconstructor.h"
#include "dualgraph.h"
#include "cut.h"

using namespace PBVP;

DCTConstructor::DCTConstructor(PBVP::Reg3DVolumeData *pData, PBVP::AugmentedContourTree *pCTree)
{
	m_pData = pData;
	m_pVolCritical = new Volume3DCritical(m_pData);
	m_pACTree = pCTree;
}

DCTConstructor::~DCTConstructor()
{}

void DCTConstructor::computeNodeAttributes(PBVP::DualNode *pDCTNode)
{
	// printf("computing DCT node attributes ....\n");
	// geometry attribute
	GeoAttribute *geoAttrib = new GeoAttribute();

	// shape attribute
	ShapeAttribute *shapeAttrib = new ShapeAttribute();

	// potential attributes
	int i;
	std::vector<PotentialAttribute *> vPotAttrib;
	for (i = 0; i < m_pVolCritical->getNumofProperties(); i++)
	{
		PotentialAttribute* potAttrib = new PotentialAttribute;
		vPotAttrib.push_back(potAttrib);
	}

	float range[2];
	range[0] = pDCTNode->fmin;
	range[1] = pDCTNode->fmax;

	// total number of tetrahedron = 6 * number of cells
	int numTetra = 6 * m_pData->numOfCells();

	bool *touched = new bool[numTetra];
	memset(touched, 0, numTetra);

	// The queue of tetrahdron to be marched
	std::deque<int> tetraQueue;
	// add the first tetra
	touched[pDCTNode->startTetra] = true;
	tetraQueue.push_back(pDCTNode->startTetra);

	float face_rang[2];

	while (!tetraQueue.empty())
	{
		int tid = tetraQueue.front();
		tetraQueue.pop_front();

		float vals[4], x0[3], x1[3], x2[3], x3[3];
		m_pVolCritical->getTetraData(tid, vals, x0, x1, x2, x3);

		Tetra *tetra = new Tetra(x0, x1, x2, x3, vals[0], vals[1], vals[2], vals[3]);
		// compute volume and surface areas
		geoAttrib->lowerArea += (float) tetra->isosurfAera(range[0]);
		geoAttrib->upperArea += (float) tetra->isosurfAera(range[1]);
		//geoAttrib->volume += (float) (tetra->innerVolume(range[1]) - tetra->innerVolume(range[0]));
		geoAttrib->volume += (float) tetra->funcIntegralDiff(1.0f, 1.0f, 1.0f, 1.0f, range[1], range[0]);
		// compute shape integral
		// nine functions: x, y, z, x^2, y^2, z^2, xy, xz, and yz
		shapeAttrib->x_Integral += (float)(tetra->funcIntegralDiff(x0[0], x1[0], x2[0], x3[0], range[1], range[0]));
		shapeAttrib->y_Integral += (float)(tetra->funcIntegralDiff(x0[1], x1[1], x2[1], x3[1], range[1], range[0]));
		shapeAttrib->z_Integral += (float)(tetra->funcIntegralDiff(x0[2], x1[2], x2[2], x3[2], range[1], range[0]));
		shapeAttrib->x2_Integral += (float)tetra->funcIntegralDiff(x0[0]*x0[0], x1[0]*x1[0], x2[0]*x2[0], x3[0]*x3[0], range[1], range[0]);
		shapeAttrib->y2_Integral += (float)tetra->funcIntegralDiff(x0[1]*x0[1], x1[1]*x1[1], x2[1]*x2[1], x3[1]*x3[1], range[1], range[0]);
		shapeAttrib->z2_Integral += (float)tetra->funcIntegralDiff(x0[2]*x0[2], x1[2]*x1[2], x2[2]*x2[2], x3[2]*x3[2], range[1], range[0]);
		shapeAttrib->xy_Integral += (float)tetra->funcIntegralDiff(x0[0]*x0[1], x1[0]*x1[1], x2[0]*x2[1], x3[0]*x3[1], range[1], range[0]);
		shapeAttrib->xz_Integral += (float)tetra->funcIntegralDiff(x0[0]*x0[2], x1[0]*x1[2], x2[0]*x2[2], x3[0]*x3[2], range[1], range[0]);
		shapeAttrib->yz_Integral += (float)tetra->funcIntegralDiff(x0[1]*x0[2], x1[1]*x1[2], x2[1]*x2[2], x3[1]*x3[2], range[1], range[0]);

		// compute the potential integrals
		for (i = 0; i < m_pVolCritical->getNumofProperties(); i++)
		{
			PotentialAttribute* potAttrib = vPotAttrib[i];
			// potential values
			float pot[4];
			m_pVolCritical->getTetraProperty(tid, pot, i);
			potAttrib->p_Integral += (float)tetra->funcIntegralDiff(pot[0], pot[1], pot[2], pot[3], range[1], range[0]);
			potAttrib->px_Integral += (float)tetra->funcIntegralDiff(pot[0]*x0[0], pot[1]*x1[0], pot[2]*x2[0],
				pot[3]*x3[0], range[1], range[0]);
			potAttrib->py_Integral += (float)tetra->funcIntegralDiff(pot[0]*x0[1], pot[1]*x1[1], pot[2]*x2[1],
				pot[3]*x3[1], range[1], range[0]);
			potAttrib->pz_Integral += (float)tetra->funcIntegralDiff(pot[0]*x0[2], pot[1]*x1[2], pot[2]*x2[2],
				pot[3]*x3[2], range[1], range[0]);
			potAttrib->px2_Integral += (float)tetra->funcIntegralDiff(pot[0]*x0[0]*x0[0], pot[1]*x1[0]*x1[0], pot[2]*x2[0]*x2[0],
				pot[3]*x3[0]*x3[0], range[1], range[0]);
			potAttrib->py2_Integral += (float)tetra->funcIntegralDiff(pot[0]*x0[1]*x0[1], pot[1]*x1[1]*x1[1], pot[2]*x2[1]*x2[1],
				pot[3]*x3[1]*x3[1], range[1], range[0]);
			potAttrib->pz2_Integral += (float)tetra->funcIntegralDiff(pot[0]*x0[2]*x0[2], pot[1]*x1[2]*x1[2], pot[2]*x2[2]*x2[2],
				pot[3]*x3[2]*x3[2], range[1], range[0]);
			potAttrib->pxy_Integral += (float)tetra->funcIntegralDiff(pot[0]*x0[0]*x0[1], pot[1]*x1[0]*x1[1], pot[2]*x2[0]*x2[1],
				pot[3]*x3[0]*x3[1], range[1], range[0]);
			potAttrib->pxz_Integral += (float)tetra->funcIntegralDiff(pot[0]*x0[0]*x0[2], pot[1]*x1[0]*x1[2], pot[2]*x2[0]*x2[2],
				pot[3]*x3[0]*x3[2], range[1], range[0]);
			potAttrib->pyz_Integral += (float)tetra->funcIntegralDiff(pot[0]*x0[1]*x0[2], pot[1]*x1[1]*x1[2], pot[2]*x2[1]*x2[2],
				pot[3]*x3[1]*x3[2], range[1], range[0]);
		}
		delete tetra;

		// consider four face-neighbor tetrahedra
		for (int i = 0; i < 4; i++)
		{
			face_rang[0] = std::min(vals[(i+1)%4], std::min(vals[(i+2)%4], vals[(i+3)%4]));
			face_rang[1] = std::max(vals[(i+1)%4], std::max(vals[(i+2)%4], vals[(i+3)%4]));
			if (!m_pVolCritical->areDisjointRanges(range, face_rang))
			{
				int neighbor_tid = m_pVolCritical->tetNeighbor(tid, i);
				if (neighbor_tid >= 0 && !touched[neighbor_tid])
				{
					touched[neighbor_tid] = true;
					tetraQueue.push_back(neighbor_tid);
				}
			}
		}
	}
	shapeAttrib->volume = geoAttrib->volume;
	pDCTNode->addAttribute(geoAttrib);
	pDCTNode->addAttribute(shapeAttrib);
	// topology attribute
	TopoAttribute *topoAttrib = new TopoAttribute;
	// TODO: this needs to be changed
	topoAttrib->innerBetti = pDCTNode->bl;
	topoAttrib->outerBetti = pDCTNode->bu;
	pDCTNode->addAttribute(topoAttrib);

	for (i = 0; i < m_pVolCritical->getNumofProperties(); i++)
	{
		vPotAttrib[i]->volume = geoAttrib->volume;
		vPotAttrib[i]->x_Integral = shapeAttrib->x_Integral;
		vPotAttrib[i]->y_Integral = shapeAttrib->y_Integral;
		vPotAttrib[i]->z_Integral = shapeAttrib->z_Integral;
		pDCTNode->addAttribute(vPotAttrib[i]);
	}

	delete[] touched;
}

GeometryData* DCTConstructor::extractDCTNodeSurface(DualNode* pDCTNode, bool upperSurface)
{
	GeometryData* surface = new GeometryData();

	// isovalue of surface extraction
	float isoval;
	if (upperSurface)
	{
		isoval = pDCTNode->fmax;
	} else
	{
		isoval = pDCTNode->fmin;
	}

	// total number of tetrahedron = 6 * number of cells
	int numTetra = 6 * m_pData->numOfCells();
	bool *touched = new bool[numTetra];
	memset(touched, 0, numTetra);


	// The queue of tetrahdron to be marched
	std::deque<int> tetraQueue;
	// add the first tetra
	touched[pDCTNode->startTetra] = true;
	tetraQueue.push_back(pDCTNode->startTetra);

	// edge cache
	EdgeCache cache;

	float range[2];
	range[0] = pDCTNode->fmin;
	range[1] = pDCTNode->fmax;
	float face_rang[2];

	while (!tetraQueue.empty())
	{
		int tid = tetraQueue.front();
		tetraQueue.pop_front();

		float vals[4], x0[3], x1[3], x2[3], x3[3];
		m_pVolCritical->getTetraData(tid, vals, x0, x1, x2, x3);

		Tetra *tetra = new Tetra(x0, x1, x2, x3, vals[0], vals[1], vals[2], vals[3]);

		int vids[4];
		m_pVolCritical->getTetraVertIds(tid, vids);
		// extract isosurf triangles from the tetra and add them to the surface
		tetra->extractSurface(vids, isoval, surface, cache);

		delete tetra;
		// consider four face-neighbor tetrahedra
		for (int i = 0; i < 4; i++)
		{
			face_rang[0] = std::min(vals[(i+1)%4], std::min(vals[(i+2)%4], vals[(i+3)%4]));
			face_rang[1] = std::max(vals[(i+1)%4], std::max(vals[(i+2)%4], vals[(i+3)%4]));
			if (!m_pVolCritical->areDisjointRanges(range, face_rang))
			{
				int neighbor_tid = m_pVolCritical->tetNeighbor(tid, i);
				if (neighbor_tid >= 0 && !touched[neighbor_tid])
				{
					touched[neighbor_tid] = true;
					tetraQueue.push_back(neighbor_tid);
				}
			}
		}
	}
	delete[] touched;
	return surface;
}


DualGraph* DCTConstructor::buildDualGraph(float fmin, float fmax, int numRanges)
{
	DualGraph *pDCTree = new DualGraph(fmin, fmax, numRanges);

	const static double MIN_THRE = 1e-6; 

	int i, j;

	// vertices generated because of cut
	CutVertex **new_verts = (CutVertex **)malloc(sizeof(CutVertex *)*(numRanges+1));

	int *nverts = new int[numRanges+1];
	unsigned int dim[3];
	m_pData->getDimension(dim);
	int max_id = std::max(int(dim[0]*dim[1]*dim[2]), m_pACTree->getMaxID());
	int cur_id = max_id;
	float factor = numRanges/(fmax - fmin);
	
	std::map<int, CutVertex*> *cut_maps = new std::map<int, CutVertex*>[numRanges+1];
	pDCTree->cut_vals[0] = fmin; 
	pDCTree->cut_vals[numRanges] = fmax;
	for (i = 1; i < numRanges; i++)
	{
		pDCTree->cut_vals[i] = fmin + i*(fmax-fmin)/numRanges;
	}
	for (i = 0; i <= numRanges; i++)
	{
		new_verts[i] = cutContourTree(pDCTree->cut_vals[i], nverts[i]);
#ifdef _DEBUG
		printf("# of cuts at value %f: %d\n", pDCTree->cut_vals[i], nverts[i]);
#endif
		for (j = 0; j < nverts[i]; j++)
		{
			new_verts[i][j].id = cur_id ++;
			(cut_maps[i])[new_verts[i][j].super_arc.id] = &(new_verts[i][j]);
		}
	}

	int *levels = new int[max_id];
	dynamic_array<int> *cpnt_sets = new dynamic_array<int>[numRanges];
	// compute the level of all critical points
	std::vector<CriticalPoint>* critPoints = m_pACTree->getSortedPoints();
	std::vector<CriticalPoint>::iterator it = critPoints->begin();
	while(it != critPoints->end()) 
	{
		int cid = (*it).id;
		float val = (*it).val;
		if (val < fmin) {
			levels[cid] = -1;
		} else if (val > fmax) {
			levels[cid] = numRanges;
		} else {
			levels[cid] = std::min(numRanges-1, int((val-fmin)*numRanges/(fmax-fmin)));
			if(val <= pDCTree->cut_vals[levels[cid]]) {
				levels[cid] --;
			}
		}
		if (levels[cid] >= 0 && levels[cid] < numRanges)
		{
			(cpnt_sets[levels[cid]]).insert(cid);
		}
		++it;
	}
	
	for (i = 0; i < numRanges; i++)
	{
		buildRange(pDCTree, i, cpnt_sets[i], new_verts[i], nverts[i], new_verts[i+1], 
				   nverts[i+1], cut_maps, levels);
	}

	delete[] levels;
	// connect the neighboring dual nodes
	for (i = 1; i < numRanges; i++)
	{
		for (j = 0; j < nverts[i]; j++)
		{
			DualNode::connect(new_verts[i][j].ngb_l, new_verts[i][j].ngb_u);
		}
	}
	//printf("***********total # of nodes************** = %d\n", nnodes);
	delete[] cut_maps;
	delete[] nverts;
	delete[] cpnt_sets;
	for (i = 0; i <= numRanges; i++)
	{
		delete[] new_verts[i];
	}
	free(new_verts);

	return pDCTree;
}

//////////////////////////////////////////////////////////////////////////
// Private functions
//////////////////////////////////////////////////////////////////////////
CutVertex* DCTConstructor::cutContourTree(float x, int &nv)
{
	const double minimal = 1e-8;

	int i = 0;

	std::vector<SuperArc>* cutArcs = m_pACTree->cut(x);
	nv = cutArcs->size();
	CutVertex* cutVerts = new CutVertex[nv];
	std::vector<SuperArc>::iterator it;
	int count = 0;
	for(it = cutArcs->begin(), i = 0; it != cutArcs->end(); ++it, i++) {
		SuperArc& arc = (*it);

		int n1 = arc.v1;
		int n2 = arc.v2;
		float x1 = m_pData->getValue(n1);
		float x2 = m_pData->getValue(n2);
		if(x1 == x2) {
			assert(x1 == x && x2 == x);
			//cutVerts[count++].super_arc = arc;
		} else if (x1 == x) {
			if (x2 > x) cutVerts[count++].super_arc = arc;
		} else if (x2 == x) {
			if (x1 > x) cutVerts[count++].super_arc = arc;
		} else {
			cutVerts[count++].super_arc = arc;
		}
		//printf("v1: %d %f, v2: %d %f, cut_val: %f\n", n1, x1, n2, x2, x);
	}
	nv = count;
	delete cutArcs;

	return cutVerts;
}

void DCTConstructor::buildRange(DualGraph *pDCTree, int ir, dynamic_array<int>& crit_pnts, 
								CutVertex* low_verts, int nl, CutVertex* up_verts, int nu, 
								std::map<int, CutVertex*> *cut_maps, int levels[])
{
	printf("building range %d\n", ir);

	int j, k;
	int total = crit_pnts.size() + nl + nu;
	DisjointSet *djs = new DisjointSet(total);
	std::map<int, int> refs;
	for (j = 0; j < nl; j++)
	{
		refs[low_verts[j].id] = j;
		djs->makeSet(j);
	}
	for (k = 0; k < crit_pnts.size(); k++, j++)
	{
		refs[crit_pnts[k]] = j;
		djs->makeSet(j);
	}
	for (k = 0; k < nu; k++, j++)
	{
		refs[up_verts[k].id] = j;
		djs->makeSet(j);
	}

	for (j = 0; j < crit_pnts.size(); j++)
	{
		int cid = crit_pnts[j];
		SuperNode& node = m_pACTree->getNode(cid);
		int ngb = node.degree();
		//SuperNode *snode = (actree->nodes)[cid];
		//int ngb = snode->degree();
		for (k = 0; k < ngb; k++)
		{
			int nid = node.getNeighborAt(k);
			if (levels[nid] == ir)
			{
				if (refs[nid] < refs[cid])
				{
					djs->unionSet(refs[cid], refs[nid]);
				}
			} else if (levels[nid] < ir)
			{	// lower neighbor crossing the sub-range
				SuperArc& arc = node.getArc(nid);
				// LL_Node<SuperArc>* sarc = snode->getArcPtr(nid);
				nid = ((cut_maps[ir])[arc.id])->id;
				assert(refs[nid] < refs[cid]);
				djs->unionSet(refs[cid], refs[nid]);
			}
		}
	}

	for (j = 0; j < nu; j++)
	{
		SuperArc& arc = up_verts[j].super_arc;
		int nid;
		if (m_pData->getValue(arc.v2) <= m_pData->getValue(arc.v1))
		{
			nid = arc.v2;
		} else
		{
			nid = arc.v1;
		}
		if (levels[nid] == ir)
		{	// >= because the critital point nid can be on the cutting edge
			djs->unionSet(refs[up_verts[j].id], refs[nid]);
		} else
		{
			// would not consider if nid is above current segment 
			if (levels[nid] > ir)
			{
				printf("ir = %d, levels[%d] = %d\n", ir, nid, levels[nid]);
				printf("critical val = %f\n", m_pData->getValue(nid));
				printf("n1 = %d val1 = %f, n2 = %d val2 = %f\n", 
					arc.v1, m_pData->getValue(arc.v1),
					arc.v2, m_pData->getValue(arc.v2));
			}
			assert(levels[nid] < ir);
			nid = ((cut_maps[ir])[arc.id])->id;
			djs->unionSet(refs[up_verts[j].id], refs[nid]);
		}
	}

	std::map<int, DualNode*> dualrefs;
	int dual_id;
	DualNode *dualnode;
	for (j = 0; j < nl; j++)
	{
		dual_id = djs->find(j);
		if (dualrefs.find(dual_id) == dualrefs.end())
		{
			dualnode = new DualNode(ir);
			dualnode->id = pDCTree->nnodes ++;
			dualnode->bl = low_verts[j].super_arc.getBettiNumbers();
			dualnode->setFuncMax(pDCTree->cut_vals[ir]);
			dualnode->setFuncMin(pDCTree->cut_vals[ir]);
			dualrefs[dual_id] = dualnode;
			pDCTree->nodes[ir].insert(dualnode);
			low_verts[j].ngb_u = dualnode;
		} else
		{
			dualnode = dualrefs[dual_id];
			dualnode->bl += low_verts[j].super_arc.getBettiNumbers();
			low_verts[j].ngb_u = dualnode;
		}
	}
	for (j = 0; j < nu; j++)
	{
		dual_id = djs->find(j + nl + crit_pnts.size());
		if (dualrefs.find(dual_id) == dualrefs.end())
		{
			dualnode = new DualNode(ir);
			dualnode->id = pDCTree->nnodes ++;
			dualnode->bu = up_verts[j].super_arc.getBettiNumbers();
			dualnode->setFuncMax(pDCTree->cut_vals[ir+1]);
			dualnode->setFuncMin(pDCTree->cut_vals[ir+1]);
			dualrefs[dual_id] = dualnode;
			pDCTree->nodes[ir].insert(dualnode);
			up_verts[j].ngb_l = dualnode;
		} else
		{
			dualnode = dualrefs[dual_id];
			///dualnode->setFuncMax((ir+1)/(float)nrang);
			dualnode->setFuncMax(pDCTree->cut_vals[ir+1]);
			dualnode->bu += up_verts[j].super_arc.getBettiNumbers();
			up_verts[j].ngb_l = dualnode;
		}
	}

	// set the function min and max of dual nodes
	for (k = 0, j = nl; k < crit_pnts.size(); k++, j++)
	{
		dual_id = djs->find(j);
		dualnode = dualrefs[dual_id];
		CriticalPoint& cpnt = m_pACTree->getNode(crit_pnts[k]);
		//CriticalPoint* cpnt = &((*(actree->p_cpmap))[crit_pnts[k]]->object);
		dualnode->setFuncMax(cpnt.val);
		dualnode->setFuncMin(cpnt.val);
		dualnode->ncrit ++;
	}

	// compute the attributes of dual nodes
	for (j = 0; j < nl; j++)
	{
		dual_id = djs->find(j);
		dualnode = dualrefs[dual_id];

		if (!dualnode->isVolumeComputed())
		{
			SuperArc& arc = low_verts[j].super_arc;
			float cut_val = pDCTree->cut_vals[ir];
			computeDCTNodeAttrib(dualnode, arc, cut_val);
		}
	}

	for (j = 0; j < nu; j++)
	{
		dual_id = djs->find(j + nl + crit_pnts.size());
		dualnode = dualrefs[dual_id];

		if (!dualnode->isVolumeComputed())
		{
			SuperArc& arc = up_verts[j].super_arc;
			float cut_val = pDCTree->cut_vals[ir+1];
			computeDCTNodeAttrib(dualnode, arc, cut_val);
		}
	}

	delete djs;
} 

void DCTConstructor::computeDCTNodeAttrib(DualNode* dualnode, SuperArc& arc, float cut_val)
{
	int v1, v2;
	int idx[3], nidx[3];
	float range[2];

	range[0] = dualnode->fmin;
	range[1] = dualnode->fmax;

	if (arc.intra_node < 0)
	{
		if (m_pData->getValue(arc.v2) < m_pData->getValue(arc.v1))
		{
			v1 = arc.v2;
			v2 = arc.v1;
		} else
		{
			v1 = arc.v1;
			v2 = arc.v2;
		}

		if (!m_pVolCritical->areNeighbors(v1, v2))
		{
			int xid;
#ifdef _DEBUG
			m_pData->vert2index(v1, idx);
			m_pData->vert2index(v2, nidx);

			printf("v1 : (%d %d %d) = %f, ", idx[0], idx[1], idx[2], m_pData->getValue(v1));
			printf("v2 : (%d %d %d) = %f, ", nidx[0], nidx[1], nidx[2], m_pData->getValue(v2));
#endif
			SuperNode& snode1 = m_pACTree->getNode(v1);
			int x1 = snode1.getNeighborTag(v2);
			SuperNode& snode2 = m_pACTree->getNode(v2);
			int x2 = snode2.getNeighborTag(v1);
			if (m_pData->getValue(x1) > m_pData->getValue(v1) && m_pVolCritical->areNeighbors(v1, x1))
			{
				xid = x1;
				v2 = xid;
			} else if (m_pData->getValue(x2) < m_pData->getValue(v2) && m_pVolCritical->areNeighbors(v2, x2))
			{
				xid = x2;
				v1 = xid;
			}
		}
		assert(m_pVolCritical->areNeighbors(v1, v2));
		// printf("v1 degree = %d, v2 degree = %d\n", (m_pACTree->nodes)[v1]->degree(),
		//   (m_pACTree->nodes[v2])->degree());
	} else {
		v1 = arc.intra_node;
		v2 = m_pVolCritical->getCutEdge(cut_val, range, v1);
		if (m_pData->getValue(v1) > m_pData->getValue(v2))
		{
			int tmp = v1;
			v1 = v2;
			v2 = tmp;
		}

#ifdef _DEBUG
		printf("v1: %d %f, v2: %d %f, cut_val: %f, range: (%f %f)\n", v1, m_pData->getValue(v1),
			v2, m_pData->getValue(v2), cut_val, range[0], range[1]);
#endif
	}
	m_pData->vert2index(v1, idx);
	m_pData->vert2index(v2, nidx);
	printf("idx : (%d %d %d) %f, nidx : (%d %d %d) %f\n", idx[0], idx[1], idx[2], m_pData->getValue(v1),
		nidx[0], nidx[1], nidx[2], m_pData->getValue(v2));
	int tid = m_pVolCritical->getTetraIndex(idx, nidx);

	dualnode->startTetra = tid;

	// dualnode->moments = data->calcMoments(range, tid);
	computeNodeAttributes(dualnode);
	//dualnode->computeAttributes(data, tid);
	// dualnode->print();
	//dualnode->volume = dualnode->moments.vol / data->getVolume();
	//dualnode->fint = dualnode->moments.Fint / data->getVolume();
	dualnode->rawVolume = dualnode->getVolume();
	if(m_pVolCritical->getNumofProperties() > 0) {
		dualnode->fint = dualnode->getFuncIntegral(); // dualnode->moments.Fint;
	}
	dualnode->setVolumeComputed();
}